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Abstract 

Swollen stacks of finite-size disc-like Laponite clay platelets are investigated within a Wigner-Seitz cell model. Each 
cell is a cylinder containing a coaxial platelet at its centre, together with an overall charge-neutral distribution of 
microscopic co and counterions, within a primitive model description. The non-linear Poisson-Boltzmann (PB) equa- 
tion for the electrostatic potential profile is solved numerically within a highly efficient Green's function formulation. 
Previous predictions of linearised Poisson-Boltzmann (LPB) theory are confirmed at a qualitative level, but large 
quantitative differences between PB and LPB theories are found at physically relevant values of the charge carried 
by the platelets. A hybrid theory treating edge effect at the linearised level yields good potential profiles. The force 
between two coaxial platelets, calculated within PB theory, is an order of magnitude smaller than predicted by LPB 
Q ' theory 

<N 

Q\ ' I. INTRODUCTION 



C3 



Electrostatic interactions between suspended mesoscopic particles or polyions play a key role in determining the 
stability, mesostructure and phase behaviour of colloidal dispersions and polyelectrolytes. Polyions in aqueous disper- 
sions may be rigid or flexible, and range in shape from spherical to rod-like or lamellar. Widely studied examples of 
*^ ■ polyions include polystyrene balls, elongated TMV (Tobacco Mosaic Virus) particles, stiff polyelectrolyte chains like 
DNA, flexible membranes and silicate clay platelets. The highly charged polyions strongly attract or repel microscopic 
counterions and coions (microions), leading to the formation of electric double layers characterised by highly inhomo- 
^ ' geneous charge distributions around the polyions. Theoretical investigations of the structure of such double layers, 
, and of their mutual effective interactions, have been mostly restricted to the simplest polyion topologies, including 
uniformly charged infinite planes or spheres. Starting with the pioneering work of Gouy jlj and Chapman on 
infinite planar double-layers, theoretical calculations are almost invariably based on the mean-field Poisson-Boltzmann 
(PB) approximation, which neglects excluded volume and Coulomb correlations between the microions, although some 
recent attempts have been made to include such correlations within a density functional (DFT) formulation, which 
may also account for discrete solvent effects ||. The non-linear PB equation for the local electrostatic potential may 
be solved analytically for a single planar double-layer, and two interacting double-layers in the salt-free case (i.e. 
in the absence of coions) jl],pL0; in the presence of salt, the one-dimensional problem of interacting charged planes 
is easily solved numerically HT However, PB theory becomes increasingly difficult to handle for more complicated 
geometries, and must generally be linearised to become tractable. A well-known application of linearised Poisson- 
Boltzmann (LPB) theory is the calculation of the effective interaction between colloidal spheres carrying a uniform 
surface charge, leading to the screened Coulomb potential of Derjaguin, Landau, Verwey and Overbeek (DLVO). 
LPB theory is, however, strictly speaking valid only provided the local potential energy felt by the microions is ev- 
erywhere small compared to the thermal energy k B T. This condition is rarely met in the immediate vicinity of the 
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highly-charged polyions, where the Coulomb energy becomes large compared to k B T, so that non-linearities become 
crucial and full PB theory should be used to determine the concentration profiles of the microions. The deficiencies of 
LPB theory are frequently patched up by introducing the rather vague concept of counterion condensation, leading to 
a Helmholtz-Stern layer of counterions highly bound to the polyions and henceforth to an effective valence Z e ff 
of the latter, significantly reduced in magnitude compared to the nominal valence Z. It should be stressed that coun- 
terion condensation is a well-defined concept only within the cylindrical geometry of an infinite, uniformly charged 
thin rod |9|, |l0| . In all other geometries, effective valences Z e ff can only be defined within some phenomenological 
convention |llj|l2| . In the case of spherical charge-stabilised colloids some recent direct measurements of the effective 
pair potential between polyions in the bulk of a suspension point to the validity of the functional form and range of the 
DLVO potential, provided Z e ff is considered to be an adjustable parameter, varying within a physically reasonable 
range 

This paper deals with electric double- layers in lamellar stacks of uniformly charged finite platelets, considered as a 
model for swollen clays. Due to the finite size of the platelets, edge effects come into play so that the problem ceases 
to be one-dimensional, as would be the case for stacks of infinite planes. For simplicity, the clay platelets are assumed 
to be disc-shaped and coaxial, while they are equally spaced within an infinite stack. Fairly monodisperse disc-shaped 
platelets are synthetised as Laponite |l4[ , a model clay which has received much recent experimental and theoretical 
attention. In particular a Wigner-Seitz (WS) model for such stacks has been thoroughly investigated within LPB 
theory p5| , ^6[ . Under most physical circumstances, the conditions of validity of LPB theory are far from being met in 
swollen clays, and consequently the much more difficult problem of non-linear PB theory for stacks of coaxial charged 
discs is being addressed in this paper. 

The key innovation of the present work is that the solution of the non-linear partial differential equation for given 
boundary conditions is not sought within a standard finite difference scheme [ p"7| , but rather by expressing the solution 
in terms of the appropriate Green's function. The latter is obtained analytically for the present cylindrical geometry, 
in the form of a Bessel-Dini series, and the resulting non-linear integral equation is solved numerically by a very stable 
iterative procedure. 

The remainder of the paper is organised as follows. The Wigner-Seitz model, appropriate for stacks of coaxial 
charged discs, is defined in section II. The basic assumptions and equations of PB theory are laid out in section 
III, and the Green's function methodology for solving the non-linear PB equation with cylindrical WS boundary 
conditions is presented in section IV. A hybrid PB/Debye-Hiickel theory, whereby edge effects are treated to linear 
order in the deviation of the potential from its counterpart for infinite planar geometry, is formulated in section V. 
Numerical results for potential and concentration profiles, and for the resulting osmotic properties, are presented in 
section VI. The force acting between two finite, coaxial platelets is calculated in section VII, while conclusions are 
drawn in section VIII. A preliminary account of parts of this work was published elsewhere (lSy. 



II. WIGNER-SEITZ MODEL FOR SWOLLEN CLAYS 



Consider a stack of coaxial, infinitely thin, disc-shaped clay platelets, with an average spacing H = 2h. Each 
platelet has a radius ro, and carries a uniform charge density a = —Ze/irrQ, where Ze is the total charge on a 
platelet. For Laponite, ro w I50Aand Z « 1000 Stacks are assumed to fill space in a columnar arrangement, 
with the normals to the platelets all pointing in the same direction. Each platelet is then placed, on average, at 
the centre of a Wigner-Seitz cell, of volume v = V/N, where V is the total volume of the sample and N the total 
number of platelets. In a hexagonal columnar array, the topology of the WS cell would be a prism of height H and 
hexagonal basis parallel to the platelet. The WS model thus describes a regular three-dimensional array of platelets 
with perfectly aligned axes, and is thus appropriate to describe moderately swollen stacked clays, but would not be 
appropriate for dilute dispersions of clay particles. In view of the circular shape of the platelet it is reasonable to 
replace the WS prism by a cylinder of identical volume and height. In fact it was shown explicitly in Ref. ]l5| that, 
at least within LPB theory, the results are very insensitive to the precise shape of the WS cell. Compared to previous 
studies of stacked clay platelets, assumed to be infinite planes |5| , the present study takes into account edge (or finite 
platelet size) effects. 

The WS model for a swollen clay which will be examined in this paper reduces hence to a single charged circular 
platelet placed at the centre of a coaxial cylinder, together with N + and A_ monovalent counterions and coions. 
Overall charge neutrality requires that 

N + = A_ + Z. (1) 

The molecular nature of the solvent (water) inside the cell will be ignored, i.e. the solvent is treated as a continuum 
of dielectric constant c (primitive model of electrolytes). The height H and radius 7? > /'o of the cylinder must be 
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such that ttR 2 H = v, where the volume v is determined by the macroscopic concentration of clay platelets. The latter 
does not determine the aspect ratio H/ R, which will be determined by the condition that it minimizes the total free 
energy of the microion distribution within the cell. The task is now to determine the total electrostatic potential (p(r) 
throughout the cell, and to derive from it the concentration profiles p + (r) and p _ (r) of counterions and coions, and 
the resulting osmotic properties, including the free energy, for a given volume v and aspect ratio H/R. In keeping 
with the physical meaning of the WS cell as representing the cage formed by the neighboring platelets, the component 
of the electric field E normal to the surface £ of the cylinder will be assumed to vanish at each point of that surface, 
i.e. the following Neumann boundary condition will be imposed throughout 

[n(r)-V] reE¥ >(r) = 0, (2) 

where n(r) is the normal to the surface E at r. The task thus defined will be carried out within PB theory in the 
following sections. 



III. POISSON-BOLTZMANN THEORY. 

The electrostatic potential throughout the WS cell satisfies Poisson's equation 

9 , > 47T , . , . 

VVr) = -— Pc(r), (3) 

where p c ( r ) is the total charge density, which includes the contribution of the platelet, qp(r), and that of the (mono- 
valent) microions 

p c (r)= gp (r)+e[p+(r)-p-(r)]. (4) 

In the mean-field approximation, the positions of microions are uncorrelated, so that their local concentrations are 
simply related to the local potential by the following Boltzmann distributions 

/9 ± (r)=p±e^W, (5) 

where /3 = l/k B T. It is worth remembering that Eq. ([|) is just a statement of the Euler-Lagrange equation resulting 
from the minimisation of the following free energy functional 

? [p + ( r ) ) P~ ( r )] = Fideal + ^Coulomb 

T ldeal = k B T f p a (r) [ln(A^(r)) - l] dv 

/T/ J JV 
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^Coulomb = ^ J Pc{r)(p(r)dr, 



(6) 



keeping in mind that p c (r) and ip(r) are linearly related by Poisson's equation. 

The prefactors p^ in Eq. (|^) have no physical significance when considered separately. However, their product 
Pq Pq is uniquely determined by the imposed physical conditions. In a canonical description, N + and are fixed by 
the macroscopic concentrations n + = N + /v and n_ = N-/v of counterions and coions; the latter coincides with the 
concentration n$ of added salt, i.e. = ng. The electroneutrality condition ([!]) then fixes n + , and the prefactors 



Pq are determined by the normalisation constraints 



n± = — [ p ± (r)dr. 

v Jv 



(7) 



Alternatively, if the suspension is in osmotic equilibrium with an ionic solution of concentration n' s acting as a 
reservoir, a semi-grand-canonical description is in order. In this case the chemical potential of co and counterions 
must be identical in the suspension and in the reservoir, leading to the condition @] 

PoPo = (n' s f . (8) 

A key quantity describing the solution is then the Donnan ratio n s /n' s of salt concentrations in the cell and in the 
reservoir. Substitution of Eq. (||) into Eqs. (Q) and (||) leads to the closed non-linear partial differential equation for 
the potential <p(r) 
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W(r) = ~ (q P (r) 



Po 



(9) 



This PB equation must be solved subject to the Neumann boundary condition (|2|) on the surface of the WS cylinder. 
The linearised version of Eg. (m , where the Boltzmann factors are expanded to first order in ip (LPB theory) was 
solved analytically in Refs. []L5j,[16|, under the same boundary conditions. The non-linear problem is re-formulated in 
terms of the appropriate Green's function in the next section. 



IV. GREEN'S FUNCTION METHODOLOGY 

The infinite dilution PB problem of an isolated platelet immersed with countcrions in an electrolyte solution has 
been solved numerically using a finite difference scheme [ fl7| |. However, this method is not well adapted to the case 
of a charged platelet confined to a WS cell. Here a semi-analytic method is proposed, whereby the PB equation is 
transformed into a non-linear integral equation, which can be solved iteratively. This approach requires the knowledge 
of a Green's function G(r, r') such that its 11 convolution" product with the charge density p c (r), defined by Eqs. (||) 
and (0), namely 
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(G*p c )(r) 



47T 

e 



f G(r,r') Pc (r')dr' 



(10) 



satisfies Poisson's equation (||) and the Neumann boundary condition (g); (G * p c ) coincides then with the required 
potential (p(r). Note that the integral in Eq. (|l0| ) is not, strictly speaking, a convolution since G is not a function of 
the relative position r r' only. 

By definition, the Green's function is a solution of the following linear Poisson's equation 



V 2 r G(r,r')=5(r-r'), 



(11) 



Integrating both sides of Eq. (|Tl|), and using Gauss' theorem, 

/ V 2 T G(r,r')dr = [ n(r) • V r G(r, r') dS 
= I 5(r-r')dr= 1, 

J V 



it becomes clear that the Green's function satisfying ( |11|) cannot obey the same boundary condition 
potential, but rather 



[n(r)-VU E G(r,r') 



1 



5(E) 



(12) 
as the 

(13) 



where 5(E) is the total area of the WS cell. The boundary condition (|T3|) precludes the possibility of expanding G in 
a Bessel-Dini series jl9]], similar to that used for the solution of the LPB problem in the same cylindrical geometry 
p5|,E6|. The difficulty may be overcome by adopting one of three possible routes. 



0. 



a) Eq. ( |1 1| ) is solved subject to the boundary condition 

[n(r)-V] reS (G* Pc ) (r) 

An explicit example is given in Appendix A. 

b) A modified Green's function G B is sought, which satisfies the following Poisson's equation 

Vl G B (r,r') = S(r-r') + B{r), 

where B(r) is an arbitrary neutralising background charge distribution, such that 

Bdr = -1. 



(14) 



(15) 



(16) 
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G s (r,r') can now be required to satisfy the same boundary condition (||) as the potential. It is easily verified 
that —4tt(G b * Pc)/t indeed satisfies Eq. (||), together with the appropriate Neumann boundary condition (0), 
since 





4tt r 






e Jv 



-% c (r)--S(r) f p c {v')dv' 



-in 



Pc(r) 



(17) 



where overall charge neutrality of the WS cell has been used; similarly 

[n(r)-V] reS (G B * Pc )(r) = / { [n(r)- V] r£S G e (r, r')} Pc {v') dr> = 0. 



(18) 



in view of the boundary condition imposed on G . 
c) The bare Laplace operator in Eq. ( jTl"| ) is replaced by a dressed or screened operator, and the solution G K of 

(V? - K 2 )G(r,r') =<5(r-0 (19) 

can then be made to satisfy the Neumann boundary condition (||) for any non-zero value of the inverse length 
k. The solution of this linear problem may be obtained analytically along the same lines as those leading to 
the potential within LPB theory |]L^Jl6| , as sketched in Appendix A. Subtracting K 2 cp(r) from both sides of the 
non-linear PB equation (|9|), it is then straightforward to check that the solution to that equation, subject to the 
proper boundary condition (H), satisfies the non- linear integral equation 



tp(r) 



4-7T 
6 



G 



K (r,r') (gp(r') + . 



p + c -(3e V (r') _ p - e /3ep(r') 



+ K i± ip{r >))dr>. 



(20) 



Since G K (r,r') is known analytically, Eq. ( |20|) may be solved numerically by an iterative procedure, starting 
from an initial guess for ip(r). 



Details for the analytic solutions of Eqs. (g_lj), (|15|) and fll9|), subject to the appropriate boundary conditions, are 
given in Appendix A. In practice the iterative procedure was implemented numerically using route c). Using the 
results from Appendix A, the required Green's function may be written in the form of the following Bessel-Dini series 



G K (r,r') = ^C±(0,r')Jo fy n ^) cosh 



n>l 



Z 



An. 



(21) 



where r = (r,4>,z) (cylindrical coordinates), the signs + and — correspond to the situations z > z' and z < z' 
respectively, y n is the n th root of Ji(y) — 0, Jo and J\ are the zeroth and first order cylindrical Bessel functions, and 
Ki 2 = (vl + k 2 R 2 )/R 2 . The coefficients C^(0,r') are given by 



^2A 



8{<t> -</>') Jo(y n r'/R) 
n R 2 smh[2h/A n ] J 2 {y n ) 



cosh 



h±z' 



(22) 



For an infinitely thin platelet of radius ro, the charge density qp(r), which provide the source term in Poisson's 
equation, may be simply expressed as 



q P (r) =aO{r - r) 5(z), 



(23) 



where O is the Heaviside function. The corresponding contribution to (p(r) in Eq. ( po|) may then be evaluated 
analytically 



r. OO 

/ dT'G*(ry)q P (r') = -<j^Y, A " 



J\ {Vn ra/R) Jo (y n r/R) 



y n sinh [h/A n ] Jo(y n 



cosh 



A n 



(24) 
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The non-linear integral equation (GO) was then solved for </j(r) using an iterative Picard method with underrelaxation 



At a particular iteration step i, the result from the previous iteration, </3 l_1 (r), and the analytic result (24) are 
used as input to compute the r.h.s of Eq. (p0|). The resulting potential ip(r) is then used to produce the input for the 
next iteration step by mixing it with Lp l ~ l {v) according to 

<p* (r) = a (r) + (1 - a) ip{v) (25) 

where < a < 1; typically a w 0.9. The iterative cycle is repeated until the relative difference \ip l — y> 8_1 |/|(/J l | at the 
centre of the cylindrical WS cell (which coincides with the centre of the disc) becomes smaller than a pre-set value, 
chosen to be 1 -5 in practice. 

1 /2 

The arbitrary inverse length k is chosen in the range [k d /10, 10 k d ], where k d = (87rng e 2 /(ek B) T) 1 is the inverse 
Debye length in the reservoir, in the case of calculations carried out in the semi-grand-canonical ensemble. The 
prefactors and Pq in Eq. (Eh are fixed during the calculation at values obeying the constraint (||). Note that 
the electroneutrality condition (nl) is not obeyed until the iterations have converged to the actual solution y(r), thus 
providing an additional global convergence test. If the calculations are carried out in the canonical ensemble, n$ = n_ 
is fixed, and k d = [87rnse 2 /(e/c B T)] 1 / 2 . The concentration of counterions is fixed by the electroneutrality condition 
(|l|) and the prefactors p$ and p§ are determined by Eq. (Q) at each step of the iteration, so that electroneutrality 
holds throughout the convergence process. 

The converged solution <p(r) must be independent of the particular choice of the auxiliary variable k. Cylindrical 
symmetry implies ip(r) — tp(r,z), which is calculated on a 2-dimensional grid n r x n z spanning half the cylinder. 
In practice n r — 240 grid points were used to cover the interval [0, R], which provides sufficient resolution for the 
representation of the Bessel functions, and n z = 100 points for the interval [0,h]. 60 terms were retained in the 
Bessel-Dini series (^)- Starting from the initial guess ip(r) = 0, convergence is generally achieved in about 50 to 200 
iterations. The number of iterations needed to achieve a relative accuracy of 10 -5 can be drastically reduced if the 
electrostatic potential obtained from the hybrid theory, introduced in the next section, is used as input. 



V. A HYBRID POISSON-BOLTZMANN/DEBYE-HUCKEL THEORY 

The hybrid Poisson-Boltzmann/Debye-Hiickel (PB/DH) theory was developed for the problem of a clay platelet in 
a cylindrical WS cell as a first attempt to go beyond LPB p5| . Within this symmetry the potential ip(r) — (p(r,z) 
can be expanded in a Bessel-Dini series [ p^[ 

( 3ertz,r) = Y / Mz)Jo( M ^), (26) 

n— 1 

which factorizes the dependence of the potential on the distance z along the axis of the cylinder, and the radial 
distance r = (x 2 + y 2 ) 1 / 2 . The hybrid PB/DH approach is an attempt to treat edge effects in a perturbative way, 
whilst keeping the non-linear PB description in the limit of infinite platelets. 

Substituting ( |26| ) into Eq. (^|), the leading term Ai(z) in the expansion is exponentiated, but the exponential is 
linearised with respect to the remainder of the series (n > 2) 

p ± = p±e™ {z) (lTY,M*)Jo (Vn^fj ■ (27) 

The first term in this expansion, n = 1, corresponds to the solution of the PB equation for an infinite charged plane in 
a WS slab of height 2h, with an effective surface charge a' — a(ro/R) 2 . The terms of order n > 2 may be considered 
as providing an estimate of the correction to the infinite charged plane limit (ro — > R) due to the finite size of the 
platelets J15[ . 

The concentrations of co and counterions are determined by 



nr = po' 



- ± \ [ e T A ^dz. (28) 



2h 



h 



It is therefore impossible to impose p^ = n + and p$ = n~ simultaneously. As a consequence, in the presence of added 
salt, calculations are preferably performed in the semi-grand canonical ensemble, at fixed reservoir salt concentration. 
This can be achieved with the choice 
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Po = Po = n 's- 



(29) 



For z ^ the semi-linearised Poisson-Boltzmann equation becomes 

~ . . 47re 
V>(r) = - — x 

OO j 

PU- M(Z) ~P*z M(z) - (pte- A ^+Po e Al «) J2M*)Jo ■ (30) 

n>2 J 

The boundary condition associated with the presence of the uniformly charged platelet at z = 0, r < tq is 

( d<p(z, r) \ 2na ^( \ f oi\ 

The electric field normal to the top and bottom of the cylinder vanishes, i.e. 

/ dip{z, r 



\ dz 



= 0. (32) 

z=±h 



Using ( p9| ) in ( p0| ) and then projecting on the basis of the zeroth order Bessel functions J leads to the following set 
of differential equations 

d 2 A 1 (z) , . 

— — — = bttIbti a smh A\ (z) 
dz 1 

d2 ^ Z) - {^) 2 A n {z) = 8nl B n' s A n (z) smh A^z) ,n>2 (33) 

with the boundary conditions at z = 

' dA-\_(z)\ 1 /Tq 



/ z=0+ b \ R 



dA n (z)\ 2 /r \ 1 /y n r 



and at z = ±/i 



dz / Z=Q+ b\RJ y n Jl{y n ) 
( dA n (z) 



\ dz 



= , 7i > 1 (35) 

z=±h 



where Ib — /3e 2 /e is the Bjerrum length and b = e/2iralB the Gouy length. 

This set of differential equations can be solved analytically in the salt-free case, for a WS cylindrical cell of infinite 
height pi (h — > oo, i? finite). A numerical solution must be sought otherwise. The first equation in ( |33| ) is identical 
to the non-linear PB equation for the ion distribution in a 1-dimcnsional cell model where each platelet is assumed 
to be an infinite uniformly charged plane confined with its monovalent co and countcrions to a slab. In order to solve 
this equation we follow the prescription given in Appendix B of || . A first order differential equation can be obtained 
after integrating once, using the boundary condition at z = h 



dAi(z) 
dz 



= — J 16nlBn' s (cosh^4i(z) — cosh Ai(h)). (36) 



The equation above is solved numerically using a 4th order Runge-Kutta algorithm. Different guesses for A\ (h) lead 
to different solutions A±(z). An underestimated initial guess for Ai(h) is increased until a solution A\(z) that verifies 
the boundary condition at z = is found. 

The second equation in ( |33"| ) can also be solved by first reducing it to a first order differential equation and then 
using an appropriate numerical algorithm for its integration. Defining a function 

L(z) ee 8nl B n' s smhAi(z) - (^-) 2 , (37) 
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one obtains 



d 2 A n (z) 
dz 2 



L(z)A n (z)=0, 



(38) 



which is a stiff differential equation if L{z) is large. This can be reduced to a first order differential equation by 
applying a Ricatti transform 



which leads to 



r)(z) 



drj{z) 
dz 



1 dA n (z) 
A n (z) dz 



L(z)- V (zf 



(39) 



(40) 



The differential equation above is again solved using a 4th order Runge-Kutta algorithm, and forced to satisfy the 
boundary condition at z = h re- written as 



n(h) = o. 



Integration of (B9J) leads to 



A n (z) = Cexp 



n(z')dz' 



(41) 



(42) 



The integral above can be evaluated numerically while the integration constant C is determined by the boundary 
condition at z = 



C = 



1 



Ty(O) 



dA n (z) 
dz 



(43) 



z=0 



The Bessel-Dini coefficients A n (z) in the expansion ( f26[ ) of the electrostatic potential can be computed numerically 
using the procedure just described. This involves solving as many independent differential equations as terms kept in 
the expansion, ie one for each coefficient. The computational cost is considerably reduced compared to the solution 
of the non-linear PB equation using the iterative scheme prescribed in the previous section. 



VI. RESULTS 



A. Potential and density profiles 

Figs. |l|a-c compare the dimensionless electrostatic potential $(z,r) = f3eip(z,r) as obtained within LPB Jl5| , |l6| | , 
PB/DH and PB. The PB/DH potential was computed by first solving ( |33"| ) for the coefficients A n {z) and then using 
these in the Bessel-Dini expansion ( p6| ) for the electrostatic potential. Only the projections (z — 0,r) and (z,r = 0) 
are shown. The results are for an aqueous solution of clay discs of radius ro = 150Aat temperature T = 300K. The 
aspect ratio of the cylindrical WS cell is h/ro = 1.25 and the concentration of clay n = 5x 10~ 3 M (1M = 1 mol 
dm -3 ). Although in the present model the clay platelets are assumed to be infinitely thin, a packing fraction can be 
defined by assining a finite thickness d to the platelets; d can be chosen to be lnm, as appropriate for Laponite [ pi) ]. 

The PB and PB/DH results in Fig. |] are semi-grand canonical calculations for n' s = 5 x 10~ 3 M and Z = 5, 100 
and 1000. The Donnan ratios ns/n' s in Table [j] compare well, indicating that the concentration of added salt in the 
WS cell is approximately the same in PB and PB/DH, for a fixed Z. Unlike other thermodynamic quantities (see 
next section), this macroscopic property is well estimated within the PB/DH theory. 

The LPB results of Fig. |l| have been produced for the same salt concentrations ns in the cell as those obtained 
within the semi-grand canonical PB theory (see Table |) . Once ns is known, the Donnan ratio reads within LPB |2l]] 

ns/n' s = (1 + Zn/n s y 1/2 . (44) 

The corresponding data are collected in Table |[ The Donnan ratio in LPB is smaller that in PB and the difference 
increases with Z. For a given salt concentration in the reservoir, LPB theory underestimates the concentration of 
added salt in the WS cell. 
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Turn next to the comparison of the electrostatic potentials in Fig. H. Although at very low surface charge [Z = 5) 
all theories produce an identical electrostatic potential, as soon as Z is as large as 100, it becomes obvious that the 
linearisation of the Boltzmann factors is no longer an acceptable approximation. However, at such high surface charge 
densities, the PB/DH still provides a very good approximation for ip(z, r). At physically relevant charge densities, e.g. 
Z = 1000, the LPB approximation gives rise to an electrostatic potential substantially different from that obtained 
within PB or PB/DH. The PB/DH theory still yields a good approximation for the potential except near the surface 

of the platelet, for r ~ tq. This illustrates the breakdown of linearisation to account for finite size corrections. 

Up to a rescaling factor, the shapes of PB and LPB potentials are similar (Fig. [I]). For a given charge density 
on the platelets (fixed Z), one can then define an effective charge Z* from a maximum likclyhood criterion between 
ip PB (Z) and ip LPB \Z*)\ for Z = 1000 and the parameters of Fig. |J the corresponding effective charge is Z* = 375 
(see Fig. ||). As expected, Z* < Z, which means that the counterions recondense within a thin layer around the 
particle until the electric potential near the surface is lowered to a value of a few k B T. Fig. || shows that the resulting 
agreement is only qualitative and that the renormalized LPB potential can neither account for the edge effects nor 
for the behaviour along the z axis. However, in the vicinity of the WS surface and far from the platelet where the 
variations of the potential are weak, PB theory can be linearised and the potential retains an LPB form, provided 
the bare charge is replaced by an effective charge Z* . In the simpler spherical geometry, given the non-linear PB 
electric potential, Alexander et al. [Tl| ] proposed to determine the effective parameters for the linearised theory by 
matching to the far-field limit of the PB solution: the resulting LPB potential, matching the PB potential up to 
its third derivative, is then an effective potential for the colloid-colloid interactions. The phase diagram associated 
with the latter effective interactions was later shown to be in good agreement with experiments for charged spherical 
colloids p^| . In the case of non-spherical colloids in anisotropic cells, the condensation of micro-ions may well be 
non uniform so that the effective particle may have a different charge distribution than the original one. Preliminary 
results confirm this scenario for clay platelets: the effective charge Z* depends on the point chosen on the WS surface 
to match the LPB and PB potentials, except of course when Z is low enough. The concept of charge renormalization 
thus seems to be inadequate for clay particles due to the non uniformity of micro-ions condensation. 



B. Osmotic properties and quadrupole moment 

Once the potential cp(r) and the concentration profiles p + (r) and p~(r) have been determined within PB, PB/DH 
or LPB theories, a number of osmotic properties may be calculated, as discussed in detail in Rcfs. |l5|,[l6). The 
Helmholtz free energy is obtained by substituting the profiles into the functional (^|), with the result 

0F = f3{U P - U c ) + {N + + AT) [In (n' s A 3 ) - l] . (45) 

where A is an irrelevant length scale, and Up and Uc are the following two contributions to the internal energy 
U = U P + Uc 

U P = ~J q P (rMr)dr, U c = t J [ p +{ T ) - p -( r )] ^( r ) dr. (46) 

Expression ( (45| ) is valid within the semi-grand canonical ensemble (which was used throughout), with the choice of 
Po = Pa = n 's- The resulting grand potential is 

0tl = PF- 2N~ In (n s \ 3 ) . (47) 

For a given clay concentration, and hence for a given WS volume v, Q and all osmotic properties depend on the aspect 
ratio h/ R (or equivalently h/rg). The equilibrium topology of the columnar stacking of platelets, as represented by the 
WS cell model, is that which minimizes as a function of h/ro, for given values of v, n' s , T and e. For convenience, 
the results shown in the figures are shifted by f2o = Z \n(n' s X 3 ), ie 

/3(fi-fio) = P(U P -Uc) + N- +N+. (48) 
The osmotic pressure n can be derived from the pressure tensor Jl5| , ^6[ 

n = \ E 7 S + ^(W E , (49) 

1 k=+- 

5] 

where (..) denotes an average over the total surface £ of the WS cell. The osmotic pressure in the reservoir is 
n = 2k B Tn' s . For a cylindrical WS cell a disjoining pressure can be denned by 
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^ = \ E P fcS ' + ^|V^ S ', (50) 
p fc=+,- 

where S" C S is the surface of the bottom and top of the WS cell. The swelling arises from the osmotic pressure 
exerted by the micro-ions and a normal spacing between the platelets (equal to the height of the cylindrical cell) is 
obtained when an axial pressure 11^ is applied. 

The a priori lowest order non-vanishing multipole moment of the charge distribution in the WS cell is the zz 
component of the traceless quadrupole tensor 

= 1 ldr[Pc(r)](2z 2 -x 2 -y 2 ). (51) 

Upon integration by parts, this may be re-expressed as a surface integral 

Ql°z =-L\ i ^)V(2z 2 - x 2 - y 2 ) ■ dS, (52) 

where Poisson's equation has been used. Numerical accuracy of the solution w(r) to the PB equation can be checked 
by comparing the quadrupole as obtained from the two equations ( pT| ) and ([52]). In the plots, is normalised by 
the quadrupole moment of the disc 



Qf z SC = \j \p(r) (2z 2 -x 2 ~y 2 )dr 



Zpr 2 

= -T 1 - (53) 

Figs, ^a-c show results for the osmotic and disjoining pressures, grand potential and quadrupole moment as functions 
of the aspect ratio h/ro, for Z = 100. Results from PB and PB/DH are compared. Within the latter theory, the 
concentration profiles may be calculated from Eq. (|27j), once the coefficients A n have been obtained; alternatively, 
the profiles may be calculated by re-exponentiating the potential determined from Eq. (p6|), according to Eq. (|^). 
Fig. |^a shows that the osmotic pressure II goes through a minimum at a well-defined aspect ratio. Within PB theory, 
the disjoining pressure coincides with the osmotic pressure at that minimum. The re-exponentiated version of 
PB/DH theory leads to results in reasonable agreement with the PB data, whereas straight PB/DH theory yields 
an unphysical osmotic pressure which is lower than the reservoir pressure Ho; note that although the corresponding 
disjoining pressure agrees well with the PB data, it intersects the osmotic pressure curve away from the minimum in 
the latter, and at a much higher aspect ratio. Fig. ||b shows that the PB results for the grand potential O go through 
a fiat minimum at the same aspect ratio as the osmotic pressure. The minimum is shifted to considerably higher 
aspect ratios within PB/DH theory and its re-exponentiated form, and does not coincide with the corresponding 
minima in the osmotic pressure curves. The PB quadrupole moments calculated from Eqs. ( |5l| ) and (^2|) coincide 
and go through zero at the value of h/ro corresponding to the osmotic pressure and grand potential minima of the 
Figs||a,b. The PB/DH results are fairly close to their PB counterparts but go through zero at an aspect ratio which 
is significantly lower than the location of the minimum of the corresponding grand potential. The re-exponentiated 
PB/DH results, on the other hand, as calculated from Eqs. ( pl| ) and ( p2[ ) are inconsistent. 

Calculations at the very low platelet charge Z = 5 lead to undistinguishable results between the PB and re- 
exponentiated PB/DH theories, as one might expect from the nearly identical potential profiles in Fig. |l[ However, 
at the physically relevant charge Z — 1000, the discrepancies between the two theories become very large, as illustrated 
in Fig. ||. 

The scenario emerging from PB theory is reminiscent of the predictions of LPB theory. In the framework of the 
latter the Helmholtz free energy goes through a minimum at the same value of the aspect ratio h/ro for which 
the disjoining and osmotic pressures are equal. The quadrupole moment Q vanishes at this same h/ro- The same 
behaviour follows from the present PB results in the semi-grand canonical ensemble, showing that = and 
n = Hi at the aspect ratio that minimises the grand potential (see Figs. || and [I|) . This illustrates the robustness of 
the equivalence between 



a thermodynamic minimization 

an electrostatic criterion (Q*°* = 0) 

a mechanical equilibrium condition (n = n a 
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• the minimization of an osmotic constraint: II is minimum at the aspect ratio where II = 11^ (see Fig. |^a and 
|a. Also holds within LPB). 

The equilibrium separation between two platelets is determined by the minimisation of the grand potential (or 
Helmholtz free energy). 

Being an approximate theory for the electrostatic potential, the PB/DH is inconsistent at Z ~ 100. However, in 
spite of all the deficiencies which the PB /DH theory exhibits when trying to compute macroscopic quantities from the 
electrostatic potential, this approximation does manage to give a good first approximation for ip(r) which compares 
well with that obtained by solving the PB equation. 



VII. FORCE ACTING BETWEEN TWO PARALLEL PLATELETS 



The Green's function methodology of section [II may be extended to calculate the potential and density profiles 
around two coaxial platelets placed symmetrically inside a WS cell, within PB theory. If the two uniformly charged 
parallel discs are placed at z — ±zq from the centre of a WS cell as shown in Fig. ||, the PB equation is still given by 
@, with the source term qp(r) now including a contribution from both platelets 

q P (r) = <re(r - r) [S(z + z ) + S(z - z )} . (54) 

The boundary conditions (j^) still applies and the integral equation (^0|) for the local potential ip(r) remains unchanged, 
with the contribution ( |24[ ) from the source term now following from ( pl| ) and (|54|), with the result 



f j f s~iK ( Is lis o r A Jl{Vnr R 

/ dr G (r, r )q P r = -2a— > A„ . r r , , I 

J v R >i y n smh[h/A 



Ji(y n r Q /R) J Q (y n r/R) 



where the function T(z) is given by 



T(z), (55) 



cosh [z /A n ] cosh [{h — z)/A n ] , z < z 

T(z) = ^ cosh [z/A n ] cosh [(h - z )/AJ , — z < z < z Q (56) 

cosh[z /A„] cosh [(h + z)/A n ] , z < —z . 

An example of a PB potential profile under conditions appropriate for Laponite is shown in Fig. ^|. 

Once <p(r) and the resulting concentration profiles p a (v) are known, the local stress tensor may be evaluated at 
each point of the WS cell 

n (r) = p(r) 1 ~^V(p<S) Vf, (57) 

47T 

p(v) = k B T p Q (r) + ^|V^(r)| 2 (58) 

«=+,— 

where I denotes the unit tensor. The force acting on platelet i £ (1,2) follows by integrating the stress tensor ( |57| ) 
over both sides, i and i , of the platelet 

rF-dS,. (59) 

In the case under consideration, of two identical coaxial discs, the force Fi = — F2 is along the unit vector z of 
the cylinder axis. The electric field E = — V<p in the vicinity of platelet i may be decomposed into continuous and 
discontinuous parts 

E = Ef + E.[ c) = ±^ a n 4 + E 2 (c) . (60) 
The force Fi may then also be expressed as 

F,=a f E 4 (c) d 2 5 (i = 1,2), (61) 
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where the integration runs over the surface of platelet i. The numerical evaluation of the force from ([n]) poses 
technical difficulties associated with the removal of the discontinuity suffered by the electric field across the platelet 
which requires some care (cf. Appendix [^) . 

Returning to the definition (|5S|), the surface integral may be transformed using the mechanical equilibrium condition 



V- n (r) = 

into the following integral over the surface S, enclosing the upper half of the WS cell (cf. Fig. 

if dS. 



(62) 



(63) 



Because n- V<^ = on the outer surface £ of the WS cell, and dip(r)/dz — by symmetry on the cross-section of the 
cylinder at z = 0, the term ® V<p gives a vanishing contribution to the force Fi, which reduces to 



Fi = - / p(r)dS. 

Js 



(64) 



This may be further simplified by noting that z-n = on the lateral surface of the cylinder (parallel to the axis), so 
that finally 



p(r) dS 



z=0 



p(r) dS. 



(65) 



z — h 



The numerical consistency between results based on Eq. (|65| ) and on Eq. ( |6l| ) (cf. Appendix [bJ) has been carefully 
checked. Explicit calculations were carried out in the limit of vanishing clay concentration (n — > 0). This was 
achieved by choosing a WS cylindrical cell large enough for the local electric field to vanish before the outer surface 
of the cell is reached. For platelets with the physical characteristics of Laponite (ro = 150A and Z = 1000), and 

salt concentrations ns ~ 10 _3 M, this condition is met with a cell volume v = 4 x 10 8 A 3 (corresponding to a clay 
concentration n = 8.3 x 10~ 6 M) and an aspect ratio h/r = 1.5. A test that the WS cell was chosen sufficiently large 
is provided by checking that the calculated concentration ns of salt in the cell coincides with the preset reservoir 
concentration n' s , as expect for vanishing clay concentration. The forces calculated within PB theory are plotted 
versus the distance between platelets in Fig. [?], and compared to the analytical prediction of LPB theory |T(| 



F*(d) = (Trro)' 



, Ana' 



f°° 1 
/ J? (x) — exp 
Jo % 



ro 



X 1 + K 2 D 7g 



(66) 



For Z = 5 the PB and LPB results coincide, but in the physically relevant case Z = 1000 shown in Fig. |?], LPB 
theory overestimates the force by an order of magnitude. This finding is consistent with the overestimation of the 
absolute value of the potential around a platelet by LPB theory, as illustrated in Fig. |l|c. Note that unlike its LPB 
counterpart, the PB force does not decay exponentially with the distance between platelets. 



VIII. CONCLUSION 



Swollen stacks of monodisperse clay platelets are conveniently modelled within a Wigner-Seitz cell representation. 
Relevant electrostatic and thermodynamic properties are derived from the inhomogeneous counterion and coion con- 
centration profiles related to the local electrostatic potential by Poisson's equation. Adopting the primitive model 
point of view for the solvent and neglecting spatial correlations between macroions leads to the closed PB equation for 
the potential ip(r) which is solved subject to Neumann boundary conditions on the confining surface of the WS cell. 
Whilst the linearised LPB version of the theory yields to analytic treatment jl5|,[l6| , the present paper concerned with 
the numerical solution of the full non- linear PB problem and of a hybrid PB /DH version of the theory. Rather than 
solving the two-dimensional non-linear partial differential PB equation on a grid, it was found that a more adequate 
and stable method is to reduce the problem to a non- linear integral equation for </?(r), involving an electrostatic 
Green's function satisfying appropriate boundary conditions. Three equivalent routes are proposed, each involving 
a Green's function which may be calculated explicitly in terms of a Bessel-Dini series, similar to that used for the 
solution of the LPB problem. 

While the present results generally confirm the qualitative trends predicted by the LPB analysis [fL5j[l(| , there are 
considerable quantitative differences at physically relevant surface charge densities of the Laponite platelets. The 
main findings may be summarized as follows 
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a) LPB theory overestimates the magnitude of the local potential by typically a factor of two for a platelet of 10 3 e. 

b) Rescaling of the platelet charge to a lower effective value to force good agreement of the PB and LPB potentials 
at the centre of the platelet fails because the shapes of the potential profiles differ significantly, particularly near 
edges. 

c) The scenario of the variation of various properties with the WS aspect ratio, for a fixed value of the cell volume, 
derived from non-linear PB theory confirms the LPB predictions: the grand potential (or the Hclmholtz free 
energy) goes through a minimum at a well-defined aspect ratio h/ro; at this same aspect ratio, the osmotic 
pressure is also at its minimum, where it coincides with the monotonically decreasing disjoining pressure, while 
the quadrupole moment of the charge distribution within the WS cell vanishes. Thus the same equilibrium 
aspect ratio is selected by thermodynamic, mechanical, osmotic and electrostatic criteria. 

d) The hybrid PB/DH theory, which is linearised with respect to edge effects yields rather accurate potential 
profiles, but rather poor thermodynamic properties; the agreement with full PB theory is improved upon re- 
exponentiation of the Bessel-Dini series for the local potential. 

e) The same Green's function methodology allows an accurate calculation of the force between the two coaxial 
platelets. The force calculated within PB theory is an order of magnitude smaller than the LPB prediction, and 
decreases with the distance between the platelets in a non-exponential fashion, up to distances of the order of 
the radius ro of the platelets. 

In future work it is planned to include microion correlations and discrete solvent (hydration) effects, by generalising 
the PB functional (||), along the lines proposed by Biben et al. ||. Generalisation of the present methodology to 
flexible charged membranes would also be helpful in describing charged soap films or smectite clay particles of lateral 
dimensions larger than those of Laponite. 
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Z = 5\ Z = 100 Z=1000 



ns/n' s 



LPB 


0.975 


0.648 


0.233 


PB/DH 


0.976 


0.720 


0.564 


PB 


0.976 


0.725 


0.574 


n s /10" 3 M 


LPB and PB 


4.878 


3.623 


2.870 


PB/DH 


4.881 


3.602 


2.820 



TABLE I. Values of the Donnan ratios n a /n' s obtained within LPB, PB/DH or PB theory, for an aqueous solution of clay 
concentration n = 5x 10 _5 M. The temperature is 300K and the radius of the discs ro = 150A. The reservoir salt concentration 
in PB and PB/DH is n' s = 5 x 10 _3 M. The concentrations of added salt in the WS cell are also indicated. 
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APPENDIX A: 



For the three routes outlined in section IV, G (G B in case b or G K in case c) is expanded in a Bessel-Dini series 



G(r, r') = Y,Gn (<M,r') Jo (vn^) 



(Al) 



Note that the coefficients of any expansion 



n>l 



(A2) 



can be obtained from the inversion relation 



fa 



R 2 J 2 o{yn) Jo 



■f(r)Jo {vn-^j dr, 



so that the Dirac distribution inside the WS cell can be cast in the form 
S(r -r') = - S(r - r') 8{<p - <f)') 8{z - z') 



(A3) 



(A4) 



Consider specifically route c. Subsitution of Eqs. (A4) and (Al) into ( [1Q[ ) yields 

l2n a 2 



fGn_A^ 2J (y n r'/R) 
dz* i? 2 n R 2 J^(y n ) 



<l>')6(z-z') 



(A5) 



where A n = Dn/R + k . With the boundary condition 



= 0, 



z=±h 



(A6) 



the solution reads 



G n (<p,z,r') = C„ cosh 



h =F z 
A„ 



(A7) 



where the superscripts + and — refer to the situations z > z' and z < z' respectively. The coefficients C„ are 
conveniently obtained by writing G n in the form 



G n (<j),z,r') = g+(z,(f>)G(z - z') + g n (z, 4>) Q(z' - z). 



(A8) 



Invoking the identity 



f(z) ^8{z - z') = -f{z) 5(z - z') + f(z') ±5{z - z') 
ciz az 



the continuity conditions obeyed by G n at z = z' follow as 



(A9) 



9n( z ') =9n( z ') 

^g + n (z') - ±g;V) = — J^y Jo (vn y S(d> - A 

leading back to Eq. ( p2|) 

An example of Green's function associated with method a can be obtained by imposing 



(A10) 



dG n 

dz 



0, for n > 2 



z=±h 



(AH) 
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while the non -Ne umann character of G reflects itself in the boundary condition obeyed by G\ , that follows from the 
resolution of (A5) with n = 0. One finds 



G 1 (z,<f >) v') = ^8(cf>-<t>')\z-z'\ 

whereas the remaining terms G n , n > 2 are the same as in and (|22]), with k = 0. 

Finally, method b can be illustrated along similar lines with the following choice of the background 



(A12) 



B 



hR 2 



(A13) 



The above density is such that the expansion (|2l|), wit h re lation ( p2[ ) and k = 0, is unchanged for n > 2 while the 
first term is now slightly modified with respect to Eq. (A.12), namely 



G 1 {z^,t') = 



1 



1 



i? 2 |Z Z 1 2hR 2 



(z 2 + z> 2 ) 



(A14) 



APPENDIX B: 



From Eq. 



the force acting on platelet 1 sitting at z = +zq can be written 

/ dip(z, r) " 



where — is for 



Zq and + is for z 



pro 

-2ira / drr 
Jo 



\ dz 



2vr 
± —a 

e 



(Bl) 



This equation allows the immediate computation of the force from the 



electrostatic potential. However, it should not be used directly since the numerical differentiation of the potential for 
z = zo and < r < tq (where dtp(r)/dz is discontinuous) is inaccurate. 
Instead, apply Gauss' theorem to (Bl) to obtain 



F = a 27TCT 

e 

Q2 _ 

= —a Zira 



\ dz 



dr ; 



f d(p(z,r) 
V dz 



Zea 
Zea 



+ 



2e ' 



(B2) 



where Q\ is the total charge (excluding the platelet) in the WS cell volume limited by zo < z < h and Qi is the total 
charge in the volume of the WS cell limited by < z < z 



Qi = 2n / dz drrp c (r) 

Jz+ Jo 
f z o f R 

Q2 = 2?r / dz drrp c (r). 
Jo Jo 

Due to the imposed electroneutrality Qi and Q2 verify 

Qi + Q2 = eZ. 



(B3) 
(B4) 

(B5) 



Eq. (|B2| ) requires the computation of the derivative of the electrostatic potential at z — zq and vq < r < R, where 
dtp(r)/dz is continuous. The differentiation can be accurately performed by fitting tp(z,r) to a polynomial in z, at 
fixed r. 
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0.0 50.0 100.0 150.0 200.0 

r or z (Angstrom) 

FIG. 1. Dimensionless electrostatic potential 3>(«,r) = (3eip{z,r) at z = versus r, and at r = versus z. The profiles 
were obtained in LPB (dotted line), in PB/DH (dashed line) and in the full non-linear PB theories (solid line). (See text for 
details.) The surface charge on the clay discs is: (a) Z = 5, (b) Z = 100 and (c) Z = 1000. The concentrations of added salt 
and the Donnan ratios ns/n' s for these results are indicated in Table I. 
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Z=1000 

LPB renormalised charge Z =375 




50 1 00 1 50 

r or z (Angstrom) 



200 



FIG. 2. Comparison between the non-linear PB potential for Z = 1000 (as in Figure [j]) and the linearised LPB potential, 
now with an effective renormalised charge Z* — 375. 
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(a) 



Z=100 


H h osmotic PB/DH 


H hdisjoining PB/DH 

x — x osmotic PB/DH re-expon. 

x x disjoining PB/DH re-expon. 

G- — ©osmotic PB 
O O disjoining PB 






0.95 



0.4 0.6 0.8 1 1.2 1.4 

h/r„ 



(b) 



Z=100 




0.4 0.6 0.8 1 1.2 1.4 

h/r 

FIG. 3. Results obtained within the PB and PB/DH theories for an aqueous solution of clay discs of charge Z — 100 and 
radius ro = 150 A . The clay concentration is n — 5 x 10" 5 M and T = 300^. The o are PB results, + are results from PB/DH 
theory and x are from PB/DH with re-exponentiation. In (a) the osmotic pressure, n, is represented by dashed lines and the 
disjoining pressure, n^, by solid lines, (b) shows the grand potential while (c) shows the quadrupole moment. The quadrupole 
obtained within PB/DH with re-exponentiation is not consistent: different results are obtained from Eqs. ([si|) (solid line) and 
© (dashed line). 
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(a) 



(b) 



Z=1000 



x — x osmotic PB/DH re-expon. 

X X disjoining PB/DH re-expon. 

G- — osmotic PB 
G disjoining PB 



2.5 

O 

C 

G 2.0 
o 

o 

c 

1.0 



0.5 



0.4 




x g x o 



0.6 0.8 1.0 

h/r„ 



Z=1000 



1.2 



1.4 



0.5 



O 

o 



-0.5 



-1.0 



X — X PB/DH re-expon., Eq. (51) 

x x PB/DH re-expon., Eq. (52) 

G OPB 



-X— x- X X X X X xx _ 



0.4 0.6 0.8 1.0 1.2 

h/r n 



1.4 



FIG. 4. Results obtained within the PB and PB/DH theory with re-exponentiation for an aqueous solution of clay discs 
of charge Z = 1000 and radius r = 150A . The clay concentration is n = 5 x 10 _5 M and T = 3007^. The o are PB 
results and x are from PB/DH with re-exponentiation. In (a) the solid lines are disjoining pressures and the dashed lines are 
osmotic pressures, (b) shows the quadrupole moment. Different results for the quadrupole are obtained within PB/DH with 
re-exponentiation, from Eqs. ( |5l| ) (solid line) and ( |52| ) (dashed line). 
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FIG. 5. Side view of the Wigner-Seitz cylinder associated with the 2 platelets problem 
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FIG. 6. Dimensionless electrostatic potential $(z,r) = f3etp(z,r) at z = versus r and at r = versus z. The profiles were 
obtained by solving the PB equation for an aqueous solution of two Laponite discs (Z = 1000, ro = 150A) at a temperature 
T — 300K. The distance between the two discs is d — 0.24ro, the height of the cylindrical WS cell is h = 1.5ro and the aspect 
ratio h/R = 0.423. The concentration of monovalent added salt in the WS cell is ns = 4.64 x 10 _3 M. 
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FIG. 7. The force between two parallel platelets of Laponite as a function of the distance separating their coaxial centres 
(F a = 47r 2 r,V 2 /<0- The solid lines are results obtained from non-linear PB and the dashed lines from LPB. The O are f° r 
n s = 0.00078M, □ for n s = 0.0046M, and the A for n s = 0.0095M. 
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